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Abstract 

Simulations are presented of the photoionisation of three dense gas clouds threaded by magnetic fields, showing the 
dynamical effects of different initial magnetic field orientations and strengths. For moderate magnetic field strengths 
the initial radiation-driven implosion phase is not strongly affected by the field geometry, and the photoevaporation 
flows are also similar. Over longer timescales, the simulation with an initial field parallel to the radiation propagation 
direction (parallel field) remains basically axisymmetric, whereas in the simulation with a perpendicular initial field 
the pillar of neutral gas fragments in a direction aligned with the magnetic field. For stronger initial magnetic fields, the 
dynamics in all gas phases are affected at all evolutionary times. In a simulation with a strong initially perpendicular 
field, photoevaporated gas forms filaments of dense ionised gas as it flows away from the ionisation front along field 
lines. These filaments are potentially a useful diagnostic of magnetic field strengths in H n regions because they are 
very bright in recombination line emission. In the strong parallel field simulation the ionised gas is constrained to 
flow back towards the radiation source, shielding the dense clouds and weakening the ionisation front, eventually 
transforming it to a recombination front. 
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1. Introduction 

Massive stars are born in cold and dense molecular 
clouds with gas temperatures of T * 15 - 50 K, but very 
early in their lives such stars reach surface effective tem- 
peratures of T e g « 30 000 - 50 000 K and emit ionising 
photons at rates of 10 47 - 10 50 s _1 [1], rapidly ionising 
their environment. For Galactic heavy element abun- 
dances, the balance of photo-heating and radiative cool- 
ing gives T « 5 000 - 10 000 K in photoionised gas fl. 
This heated region is hence a very over-pressurised bub- 
ble of ionised gas around newly formed massive stars 
(the ratio of ionised gas pressure to that in surrounding 
neutral gas w 10 2 - 10 3 ), and is known as an H n region. 
The bubble expands by driving a shock outwards J^; 
when this shock is isothermal a thin, dense, and gener- 
ally unstable shell forms at the H n region boundary. 

The interstellar medium (ISM) in molecular clouds 
is generally observed to be clumpy and filamentary 
yj, |5[] so, even in the absence of instabilities, the ex- 
pansion rate of the H n region will be uneven and angle- 
dependent. Indeed, some of the most striking and eas- 
ily observed optical nebulae are produced by recombi- 
nation radiation from ionisation fronts (I-fronts) around 



pillars of neutral gas, also known as elephant trunks |@]. 
Younger H n regions appear quite regular and spherical, 
e.g. RCW 120 [7], whereas more mature H n regions 
with ages t > 1 — 2 Myr such as the Eagle Nebula 
Ed, IC 1396 |9|. and the Carina Nebula [10] con- 
tain well-developed pillars and globules. It has been un- 
clear for decades whether ionisation front instabilities 
or pre-existing environmental inhomogeneities are pri- 
marily responsible for forming these pillars and glob- 
ules [e.g. IT, T2l; both processes should be acting, but 
possibly on different length scales. 

Magnetic fields have been suggested as the driving 
process for apparently helical and rotating structure in 
some elephant trunks 11311 . suggesting that magnetic 
fields could be dynamically important in these struc- 
tures. Only one measurement of the field orientation 
in pillars has been made so far, however, for the Eagle 
Nebula by [14]. Cometary globules are related struc- 
tures; magnetic field measurements in some of these 
. _7J show similar field alignment with the head- 
tail morphology of the globules (but see [ 18] for a coun- 
terexample). 

A magnetic field provides additional pressure support 
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to the gas (which cannot be lost by radiative cooling, 
unlike thermal pressure), and it constrains gas to flow 
along field lines (for both neutrals and ions in the ideal 
magnetohydrodynamic [MHD] limit, because of effi- 
cient collisional coupling). Pioneering calculations by 
fl9ll suggested their effects would be significant in H n 
regions; analytic models of photoionised globules in- 
cluding an approximate magnetic pressure H2Qfl showed 
that it can be important to the pressure support of such 
clouds. Extra magnetic pressure support for the Eagle 
Nebula pillars was also deduced by 12 111 by analysis of 
the gas densities and temperatures inferred from obser- 
vations. Magnetostatic turbulence B22I1 was proposed as 
a possible source of long-term pressure support, as su- 
personic turbulence decays very rapidly. Alternatively, 
the pillars may not be in pressure equilibrium but could 
still be dynamically evolving structures B231 12411 . 

The structure of magnetised ID I-fronts was stud- 
ied by I25I1 . and more recently l26Tl used 2D ionising 
radiation-MHD (IR-MHD) simulations to study slab- 
symmetric I-fronts with various initial magnetic field 
configurations, finding that strong magnetic fields dra- 
matically changed the evolution of an advancing I-front. 
The photoionisation of a dense spherical clump in a 
magnetised medium was simulated in 3D by [27]. They 
found that weakly magnetised clumps with initial field 
strength B — 59 micro-Gauss (yuG) evolved similarly to 
the purely hydrodynamic case, whereas strongly mag- 
netised clumps (B — 186 fiG) evolved in a very differ- 
ent manner. The compression and structure of the neu- 
tral globule changed and the flow of ionised gas from 
the I-front was highly constrained. For a model with a 
magnetic field perpendicular to the radiation propaga- 
tion direction, the ionised gas was confined to a dense 
ribbon or filament standing off from the globule, which 
eventually shielded the globule from much of the ion- 
ising radiation. This work was followed up by rf28ll . 
who used a similar model to to study the photoion- 
isation of multiple clumps in different configurations. 
They confirmed most of the results of 12711 . and the addi- 
tional asymmetry introduced by the multiple clumps led 
to pillar fragmentation in one case. It was also shown in 
B28I1 that the dense filament of ionised gas could be ex- 
plained quite simply by analysis of the jump conditions 
for isothermal magnetised shocks. Both studies found 
that, for a clump with a weak magnetic field oriented 
perpendicular to the radiation propagation direction ini- 
tially (hereafter referred to as perpendicular field mod- 
els), the magnetic field was swept into alignment with 
the pillar/globule during its evolution. 

In this paper two photoionisation simulations with 
magnetic fields initially aligned with the radiation prop- 



agation direction (hereafter parallel field models) are 
presented and compared to the perpendicular field simu- 
lations already described in H28I1 . for medium and strong 
initial magnetic fields. In the next section the numerical 
methods and simulation setup will be reviewed briefly; 
the results of the simulations will be presented in section 
[3j the results will be discussed in the context of obser- 
vations and previous work in section|4l and conclusions 
will be presented in section[5] 

2. Methods 

2.1. Code Description 

The IR-MHD code and simulation setup are de- 
scribed in detail in 12412811 . including results of tests of 
the algorithms, so only the essential features of the code 
are reviewed here. The ideal MHD equations (using an 
ideal gas equation of state with y - 5/3) are solved 
on a uniform 3D grid using the second-order-accurate 
(in time and space), finite-volume scheme of 112911 . with 
inter-cell fluxes calculated using the conserved variables 
Riemann solver of 1 30] (implemented as in [31]). The 
mixed-GLM method 13211 is used to control errors from 
non-zero divergence of the magnetic field. 

Microphysical processes that modify the effective 
equation of state of the gas (radiative heating and 
cooling, ionisation and recombination) are included by 
operator-splitting; in the source-term update the ther- 
mal pressure (p g ) and hydrogen ion fraction (y) are 
integrated forwards in time on a point-local cell-by- 
cell basis. Photoionisation is calculated using the 'on 
the spot' approximation, where recombinations to the 
ground state of H are ignored by assuming they result in 
re-emission of an ionising photon which is absorbed lo- 
cally. This allows one to ignore scattered radiation and 
calculate photoionisation rates based only on the direct 
photon flux from point sources (see [33] for a recent dis- 
cussion of this approximation). Attenuation of photons 
is calculated using a short-characteristics raytracer [34]. 

This algorithm is not radiation-MHD in the sense that 
radiation is not included as part of the total energy den- 
sity, and radiation pressure is not considered. Rather 
it is a photon counting scheme that tracks where ionis- 
ing photons are absorbed and deposit their energy. This 
photoionisation heating is the dominant effect of ionis- 
ing radiation for the situations we are considering. 

2.2. Simulation Setup 

The essential properties of the suite of simulations are 
stated below, in particular the radiation source, the com- 
putational grid, and the initial conditions for gas and 
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Table 1: Cloud and initial magnetic field configurations, boundary 
conditions applied, and final simulation times for the R-MHD simu- 
lations. Field strengths are quoted in micro-Gauss. Plasma parameter 
f} = p g lp m (the ratio of thermal to magnetic pressure) is shown in 
columns 5 and 6: f} n refers to the neutral gas initial conditions, and /J, 
to photoionised gas at the background density and with T = 8000 K. 



magnetic field. The reader is referred to [28] for further 
details. The simulations to be described here are models 
R5, R6, R8, and R9 from table 2 in Q; their different 
initial magnetic field configurations are given in TableQ] 
For all simulations, an ionising photon source with lu- 
minosity Qo — 2 x 10 50 photons per second is placed 
at the origin; the simulation domain is x e [1.5,6.0] 
parsecs (pc, 3.086 x 10 18 cm), and (y,z) e [-1.5, 1.5] 
pc resolved with 384 x 256 2 grid cells. Three clumps 
of gas with mass M « 28 M are placed at positions (in 
pc) (2.3, 0, 0), (2.75, 0, 0.12), and (3.2, 0, -0.12), with 
peak overdensities of 500 compared to the background 
ISM which has H number density «h = 200 cm 3 . The 
clumps have Gaussian profiles of width cr = 0.09 pc, 
and the simulation is initially in pressure equilibrium 
and static. The four simulations differ only in the ini- 
tial magnetic field strength and orientation (see TableQ}. 
All four models use zero-gradient boundary conditions. 
Models R5 and R6 have medium field strengths, mag- 
netically dominated in neutral gas, but the much hotter 
photoionised gas is thermal-pressure dominated; for R8 
and R9 both phases are magnetically dominated. R5 
and R8 have perpendicular initial field configurations, 
whereas R6 and R9 have a parallel field. 

3. Results 

Observable properties of the simulations have been 
calculated by projecting the 3D data onto a 2D plane, 
using weighted integrals along the line-of-sight (LOS). 
Three observables are plotted in Figures Q] and [2] for a 
series of times: the intensity of the Ha spectral line (n = 
3 — > 2 transition of H°), with arbitrary normalisation; 
the neutral H column density, N(H°), with units cm -2 , 
and the magnetic field orientation perpendicular to the 
line of sight. 

Evaluation of N(H°) is a simple summation of neu- 
tral gas along each column of grid cells. The Ha line 
intensity is calculated including dust attenuation but ig- 
noring radiation scattered into the LOS lEill . The gas 



temperature (T < 0.7 eV) is too low to collisionally 
populate the energy levels, so all of the Ha emission is 
from recombinations to excited states and the emissivity 
scales approximately with the recombination rate. The 
projected magnetic field is calculated using a Stokes pa- 
rameter formalism, integrated along the LOS and trans- 
formed back to a magnetic field orientation J35 , 28[] ; this 
mimics the polarimetry observations of fe.g. Il4ll . 

3.1. Simulations R5 and R6 - medium field strength 

Figure Q] shows projections along the simulation's y- 
axis, so the image plane is the simulation x-z plane. R5 
is shown at left, with the initial perpendicular field fully 
in the image plane, and R6 is shown at right, with the 
initially parallel field again fully in the image plane. The 
top two panels (at t — 0) are dark because the gas is 
initially neutral; N(H°) contours show the positions of 
the dense clumps. 

The second row shows the situation at t = 100 kyr, at 
the end of the radiation-driven implosion phase of evo- 
lution for the clump closest to the star. The main fea- 
tures of the two simulations are similar at this stage, 
with a strong radial photoevaporation (ablation) flow 
from the main I-front visible as a rim of bright Ha 
emission that weakens radially outwards because of the 
decreasing gas density. These models have B > 1 in 
ionised gas (where B = p g lp m is the ratio of thermal to 
magnetic pressure), so deviations in ionised gas flow as 
a consequence of the magnetic field are relatively mi- 
nor. The shadowed 'cometary tail' behind the clumps 
is rather different, with strong compression into a sheet 
possible in R5, whereas the parallel field in R6 resists 
compression from all sides. Neutral gas between the 
clumps is more strongly compressed in R5 than R6 for 
the same reason. 

The two models remain similar at t — 200 kyr (third 
row), but by t — 400 kyr (fourth row) substantial differ- 
ences have arisen. R6 is roughly axisymmetric but R5, 
with its strong z component of B, allows vertical motion 
(in the image plane) much more readily than horizontal 
or line-of-sight motion. Neutral gas flows in R5, driven 
in response to the photoevaporation flows, are therefore 
no longer axisymmetric and the merged clump structure 
has fragmented into two cometary globules, one much 
larger than the other. 

3.2. Simulations R8 and R9 - Strong field 

Simulations R8 and R9 are the same as R5 and R6 
but with a 3x stronger field, and small off-axis magnetic 
field components so that it is not perfectly grid aligned. 
The same plots as in Figure Q] are shown for R8 and R9 
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Figure 1: Projections through Simulations R5 (left) and R6 (right) at times (from top to bottom) 0, 100, 200, 400, 500 kyr. The image x- and 
y-axes are the simulation x- and z-axes; the line of sight is the simulation y-axis. The colour scale shows Ha intensity, integrated along the line 
of sight assuming no background sources. White lines indicate projected magnetic field orientation at the midpoint of each line, and turquoise 
contours show projected neutral hydrogen density on a linear scale. The five contours are at N(H°) = (0.2, 0.4, 0.6, 0.8, 1.0) X 10 21 cm -. Positions 
are shown on the axes labels in parsecs relative to the source; not all of the simulation is shown, and the image window moves further from the 
source from image to image. 
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in Figure now for times t = 0, 100, 200, 300, and 400 
kyr. It is immediately clear that the two models have 
completely different evolution from 100 to 400 kyr. 

As discussed in 1271 1281 . the ambient magnetic pres- 
sure confines the photoevaporation flow to a much 
smaller volume for these field strengths, and after the 
termination shock the gas flow is approximately one- 
dimensional, along field lines. For R8 (left panels) this 
leads to filaments of dense ionised gas moving toward 
the z boundaries. This dense layer shields the I-front 
at the clump surface, weakening the photoevaporation 
flow and reducing the intensity of Ha emission at the 
I-front. Indeed, the dense ribbons of ionised gas are 
the brightest regions in Ha emission. By 400 kyr the 
dense ionised layer allows significant recombination of 
the ISM at x > 4 pc (the black region in Ha). Neutral 
gas is also constrained to follow field lines and so the 
pillar structure is more anvil-shaped, elongated along 
the field direction. 

In simulation R9, the photoevaporation flow pushes 
out spherically against the magnetic pressure to some 
extent, but beyond the termination shock the postshock 
gas then follows field lines back toward the radiation 
source. This gas remains in a line between the radia- 
tion source and the dense clumps, accumulating until 
no more photons reach the dense clumps. Even at 100 
kyr (second row) significant recombination behind the 
ionised gas is clearly seen at (x,z) ~ (2.75, 0.3) pc, and 
the third clump at (3.25, -0.12) remains well shielded 
and has barely been impacted by photoionisation. A 
cometary tail has not formed to any extent because the 
shadowed region is basically incompressible (the tail in 
R8 at left is actually an edge-on thin sheet, because of 
the anisotropic magnetic pressure). From 100-400 kyr 
the Ha emission from the clump surface gets progres- 
sively weaker, and most of the emission is from the 
ionised 'plug' of gas closer to the star. Analysis of the 
gas flow at t — 400 kyr shows that the clump bound- 
ary has actually become a recombination front, where 
ionised gas recombines, cools and gets denser, joining 
the clump material. The gas velocity field has basically 
become one-dimensional over the whole simulation do- 
main, very reminiscent of the 2D simulations of IT261 
% 3]. 

The early part of R9's evolution is qualitatively sim- 
ilar to the single-clump simulation S00L of B27L fig. 8], 
although the three clumps in our calculation make the 
solution much less symmetric. Although it appears 12711 
did not follow the evolution to the point where the I- 
front in R9 switched to a recombination front, it is dif- 
ficult to make a meaningful comparison based on sim- 
ulation times because the parameters of the simulations 



are quite different. 

4. Discussion 

It is important to ask whether or not analogues of 
these idealised models exist in reality, and also to what 
extent the numerical approximations made (ideal MHD, 
limited spatial resolution, and neglecting the diffuse ion- 
ising radiation field) are affecting the solution. The first 
question can be addressed by more realistic simulations 
and detailed comparison to observations, and the second 
can potentially be addressed in laboratory experiments. 
A review of the current status of high energy-densit y ex - 
periments with photoionised plasmas is given by 13611 : 
their application to modelling of pillar formation and I- 
front instabilities is discussed by 137113811 . Extension of 
these experiments to include MHD effects requires an 
understanding of the magnetic properties of the materi- 
als used, but is in principle possible. An important is- 
sue is that the limited numerical resolution of our calcu- 
lations may artificially suppress any instabilities in the 
dynamics that are seeded on small scales (see [39]), and 
laboratory experiments could certainly test this 13811 . 

In all cases where magnetic fields have been mea- 
sured in interstellar clouds, the magnetic energy density 
is found to be comparable to kinetic and gravitational 
energy [e.g.|40j,|4lj,|42t], and larger than the thermal en- 
ergy density. Field strengths such as those in R5 and 
R6 are therefore the norm, but cases such as R8 and R9 
may be somewhat extreme. Simulations of photoioni- 
sation in a turbulent cloud l35ll showed that situations 
such as R8/R9 did not arise naturally, with the caveat 
that the initial magnetic field strength in the simulation 
was set by the specific turbulence calculation they began 
with. 

The high density of molecular pillars means that 
we are in the fully collisional MHD regime but the 
ion/electron fraction is typically very low, so imperfect 
ion-neutral coupling may be an important considera- 
tion on small scales. Multi-fluid simulations of weakly 
ionised plasmas including non-ideal MHD effects 14311 
are now beginning to test the limits of applicability of 
ideal MHD on small scales in molecular clouds. 

Infrared observations are dramatically increasing our 
understanding of massive star formation regions [e.g. 

, and millimetre data can now probe gas properties 
with high sensitivity and spatial resolution J3l ■ Together 
with recombination-line data to study ionised gas, de- 
tailed studies of the interfaces between ionised and neu- 
tral gas can now be performed Examples of 
dense sheets or filaments of photoionised gas near the 
heads of pillars have been found in NGC 6357 l44ll and 
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NGC 3603 [46]. Combined with velocity information 
in the ionised gas, and polarimetry data to constrain the 
magnetic field geometry, it would be possible to deter- 
mine if these regions are magnetically dominated, and 
so to constrain the field strength. 



5. Conclusions 

Simulations have been presented of the photoionisa- 
tion of three dense clumps threaded by magnetic fields 
of different field strengths and orientations. These are 
very idealised models, designed to show the effects of 
different initial field orientations and strengths, and they 
confirm and extend previous 2D results [26], and mod- 
els of photoionisation of a single clump in 3D M27I1 . For 
moderate initial magnetic field strengths (J3 < 1 in neu- 
tral gas, /3 > 1 in ionised gas) the initial radiation-driven 
implosion phase is not strongly affected by the field ge- 
ometry, and the photoevaporation flows are also simi- 
lar. Over longer timescales, however, the parallel field 
model R6 remains basically axisymmetric whereas the 
pillar in the perpendicular field model R5 fragments at 
late times in a direction aligned with the magnetic field. 

For stronger initial fields (J3 < 1 in both neutral and 
ionised gas), the gas dynamics in all gas phases are af- 
fected at all evolutionary times. In the parallel field 
model (R9) the I-front of the dense cloud weakens and 
becomes a steady recombination front, because pho- 
toionised gas is constrained to remain in a line between 
the radiation source and the cloud (the gas flow becomes 
almost one-dimensional). Filaments of dense ionised 
gas also develop in the perpendicular field model R8, as 
photoevaporated gas flows away from the I-front along 
field lines. These filaments are potentially a useful diag- 
nostic of magnetic field strengths in H n regions because 
they are so bright in recombination line emission. The 
absence of this emission in the recombination front in 
model R9 may also provide useful observational con- 
straints. 
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